
clear all;
clc;
disp('test result with Z formula using mean and sigma until the current point time');
load dataAnalysisQ_G1000k.mat

INSAMPLE_YR = 2008;
INSAMPLE_QR = 4;
CRITERIA = 10;

inSampleLength = (INSAMPLE_YR - 1998)* 4 + INSAMPLE_QR;
outSampleLength = length(SVD)-inSampleLength;

xx =zeros(RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
avgRollingTM =zeros(length(SVD),RATE_LIST_LENGTH-1,RATE_LIST_LENGTH); 
bin=zeros(length(SVD),RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
InsampleSVD = zeros(inSampleLength,1);
InsampleMobilityNorm = zeros(inSampleLength,1);
PHat=zeros(10,inSampleLength,RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
N = zeros(inSampleLength+outSampleLength,RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
naive1 = zeros(inSampleLength, RATE_LIST_LENGTH-1, RATE_LIST_LENGTH );
naive2 = zeros(inSampleLength, RATE_LIST_LENGTH-1, RATE_LIST_LENGTH );

filename='optimization_08q4_jun13.xls';


%To calculate the rolling avg TM and rolling credit score/bins-------------

for i = 1: NUMBER_OF_YEAR
    for j = 1: Quarter_LIST_LENGTH
        x(:,:) = durationCount(i,j,:,:);       
        xx = xx + x;
        t = (i-1)*4 + j;
        N(t,:,:)= x;
        
        for ii = 1 : RATE_LIST_LENGTH-1
            for jj = 1: RATE_LIST_LENGTH
                avgRollingTM (t, ii, :) = xx(ii,:)./ sum(xx(ii,:),2);
                cdf = sum( avgRollingTM(t, ii, jj:RATE_LIST_LENGTH) );
                bin(t,ii,jj)=norminv(cdf,0,1);
            end
        end
        
    end
end

disp('finish Rolling Avg Transition Matrix and Bin');


% To calculate credit cycle index Z, real and forecast----------------------
% We here use empirical Z to calculate w.

% Take out the in-sample data

PD = defaultFreq(1:inSampleLength);
inversePD = norminv(PD);
InsampleAvgTM (:,:) = avgRollingTM (inSampleLength, :, :);
InsampleAvgTMmobilityNorm = mean(svd(InsampleAvgTM));

for i = 1: inSampleLength
    b(:,:) = stdProbQ (i, :, :);
    InsampleMobilityNorm (i,1) = mean(svd(b));
    InsampleSVD(i,1) = InsampleMobilityNorm (i,1) - InsampleAvgTMmobilityNorm; %SVD
end

%impirical Z value

Z_PD = - zscore(inversePD); 
Z_SVD =  zscore(InsampleSVD);



% To minimize the distance to get estimate of w----------------------------
% PHat_BFS = zeros (inSampleLength, RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);

R = RATE_LIST_LENGTH-1;
C = RATE_LIST_LENGTH;
P = stdProbQ;
X(:,:) = bin(inSampleLength,:,:);
% wantedZ = Z_PD; 
wantedZ = Z_SVD;
t = inSampleLength;
W = zeros(10,1);


% L1
obj1=@(w)difL1_fixedW(w,t,P,X,wantedZ,R,C); 
[W(1,1),fval,exitflag,output] = fminbnd(obj1,0,1);
% W(1,1) = foundW;
% subplot(5,2,1); fplot(obj1,[0 1],'*');
% title('Mininization Object with L1 distance',... 
%   'FontWeight','bold');

% L2
obj2=@(w)difL2_fixedW(w,t,N,X,wantedZ,R,C); 
[W(2,1) ,fval,exitflag,output] = fminbnd(obj2,0,1);
% W(2,1) = foundW;
% subplot(5,2,2); fplot(obj2,[0 1],'*');
% title('Mininization Object with L2 distance',... 
%   'FontWeight','bold');

% WAD
obj3=@(w)difWAD_fixedW(w,t,P,X,wantedZ,R,C); % weight is est. P
[W(3,1),fval,exitflag,output] = fminbnd(obj3,0,1);
% W(3,1) = foundW;
% subplot(5,2,3); fplot(obj3,[0 1],'*');
% title('Mininization Object with WAD distance',... 
%   'FontWeight','bold');

% NAD
obj4=@(w)difNAD_fixedW(w,t,P,X,wantedZ,R,C); 
[W(4,1),fval,exitflag,output] = fminbnd(obj4,0,1);
% W(4,1) = foundW;
% subplot(5,2,4); fplot(obj4,[0 1],'*');
% title('Mininization Object with NAD distance',... 
%   'FontWeight','bold');

% SVD
obj5=@(w)difSVD_fixedW(w,t,P,X,wantedZ,R,C); 
[W(5,1),fval,exitflag,output] = fminbnd(obj5,0,1);
% W(5,1) = foundW;
% % subplot(5,2,5); 
fplot(obj5,[0 1],'*');
title('Mininization Object with SVD distance',... 
  'FontWeight','bold');

% NSD distance
obj6=@(w)difNSD_fixedW(w,t,P,X,wantedZ,R,C); 
[W(6,1),fval,exitflag,output] = fminbnd(obj6,0,1);
% W(6,1) = foundW;
% subplot(5,2,6); fplot(obj6,[0 1],'*');
% title('Mininization Object with NSD distance',... 
%   'FontWeight','bold');

% BFS distance
obj7=@(w)difBFS_fixedW(w,t,N,P,X,wantedZ,R,C); 
[W(7,1) ,fval,exitflag,output] = fminbnd(obj7,0,1);
% W(7,1) = foundW;
% subplot(5,2,7); fplot(obj7,[0 1],'*');
% title('Mininization Object with BFS distance',... 
%   'FontWeight','bold');


% D3
obj8=@(w)difD3_fixedW(w,t,P,X,wantedZ,R,C); 
[W(8,1),fval,exitflag,output] = fminbnd(obj8,0,1);
% W(8,1) = foundW;
% subplot(5,2,8); fplot(obj8,[0 1],'*');
% title('Mininization Object with D3 distance',... 
%   'FontWeight','bold');

% D1
obj9=@(w)difD1_fixedW(w,t,P,X,wantedZ,R,C); 
[W(9,1),fval,exitflag,output] = fminbnd(obj9,0,1);
% W(9,1) = foundW;
% subplot(5,2,9); fplot(obj9,[0 1],'*');
% title('Mininization Object with D1 distance',... 
%   'FontWeight','bold');

% D1sqr
obj10=@(w)difD1sqr_fixedW(w,t,P,X,wantedZ,R,C); 
[W(10,1),fval,exitflag,output] = fminbnd(obj10,0,1);
% W(10,1) = foundW;
% subplot(5,2,10); fplot(obj10,[0 1],'*');
% title('Mininization Object with D1sqr distance',... 
%   'FontWeight','bold');
disp('finish w estimation');

% To calculate in-sample fitted matrix and error(Mean Absolute Error)-------------------

for i= 1: CRITERIA
    
  PHat(i,:,:,:) = calculatePHat_fixedW(W(i,1),t,X,wantedZ,R,C); % PHat = f(X,What,Z)
    
end


for i = 1: t
if i == 1
    naive1(i,:,:) = avgRollingTM (i,:,:); % take average as naive1 benchmark
else
    naive1(i,:,:) = avgRollingTM (i-1,:,:);
end
end

for i = 1: t
if i == 1
    naive2(i,:,:) = stdProbQ (i,:,:);% take previous as naive2 benchmark
else
    naive2(i,:,:) = stdProbQ (i-1,:,:);
end
end

e = zeros(CRITERIA,inSampleLength,3);
Error = zeros(CRITERIA,3);

for i = 1:CRITERIA
    
    for j = 1:t
    x1(:,:) = PHat(i,j,:,:);
    x2(:,:)= naive1(j,:,:);
    x3(:,:) = naive2(j,:,:);
    y(:,:) = stdProbQ(j,:,:);
    
    e(i,j,:) = calculateMAE(i,x1,x2,x3,y,R,C);
      
    end
    
    Error (i, 1) = mean (abs (e(i,2:inSampleLength,1)));
    Error (i, 2) = mean (abs (e(i,2:inSampleLength,2)));
    Error (i, 3) = mean (abs (e(i,2:inSampleLength,3)));
    
end
disp('finish error calc');

% visual check est. TM-----------------------------------------------------

% first, we look at cell prob. 
% criteria D3
g1 (:,1)= PHat(8,:,1,1); g1(:,2) = naive1(:,1,1); g1(:,3) = naive2(:,1,1); g1(:,4)= stdProbQ (1:length(PHat),1,1);% rank 1.5 - 1.5 
g2 (:,1)= PHat(8,:,3,3); g2(:,2) = naive1(:,3,3); g2(:,3) = naive2(:,3,3); g2(:,4)= stdProbQ (1:length(PHat),3,3);% rank 2.5 - 2.5
g3 (:,1)= PHat(8,:,6,6); g3(:,2) = naive1(:,6,6); g3(:,3) = naive2(:,6,6); g3(:,4)= stdProbQ (1:length(PHat),6,6);% rank 3 - 3
g4 (:,1)= PHat(8,:,5,8); g4(:,2) = naive1(:,5,8); g4(:,3) = naive2(:,5,8); g4(:,4)= stdProbQ (1:length(PHat),5,8);% rank 3.5 -D
g5 (:,1)= PHat(8,:,6,8); g5(:,2) = naive1(:,6,8); g5(:,3) = naive2(:,6,8); g5(:,4)= stdProbQ (1:length(PHat),6,8);% rank 4 -D
g6 (:,1)= PHat(8,:,7,8); g6(:,2) = naive1(:,7,8); g6(:,3) = naive2(:,7,8); g6(:,4)= stdProbQ (1:length(PHat),7,8);% rank 4.5 -D

% criteria SVD

g6_svd (:,1)= PHat(5,:,7,8); g6_svd(:,2) = naive1(:,7,8); g6_svd(:,3) = naive2(:,7,8); g6_svd(:,4)= stdProbQ (1:length(PHat),7,8);% rank 4.5 -D
g2_svd (:,1)= PHat(5,:,3,3); g2_svd(:,2) = naive1(:,3,3); g2_svd(:,3) = naive2(:,3,3); g2_svd(:,4)= stdProbQ (1:length(PHat),3,3);% rank 2.5 - 2.5

s =  zeros(t,4);

for i = 1:t
    x1(:,:) = PHat(5,i,:,:); x2(:,:) = naive1(i,:,:); x3(:,:) = naive2(i,:,:); y(:,:) = stdProbQ (i,:,:);
    s(i,1) = mean (svd(x1));
    s(i,2) = mean (svd(x2));
    s(i,3) = mean (svd(x3));
    s(i,4) = mean (svd(y));
end

figure
plot (s); %([g1,gg1,ggg1]);%,'-*'); plot (gg1,'-rs'); 
set(gca,'XTick',1:4:length(g1));
labels = quaterlabels(1998, length(g1));
labels = labels(1:4:length(g1));
set(gca,'XTickLabel',labels);
xlabel('Quarter');
legend('model','benchmark1-avg','benchmark2-previous','empirical','Location','NorthWest');

% then we look at metrics

  

 
